Quantum Simulations of Classical Annealing Processes 
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We describe a quantum algorithm that solves combinatorial optimization problems by quantum 
simulation of a classical simulated annealing process. Our algorithm exploits quantum walks and 
the quantum Zeno effect induced by evolution randomization. It requires order 1/vo steps to find 
an optimal solution with bounded error probability, where 5 is the minimum spectral gap of the 
stochastic matrices used in the classical annealing process. This is a quadratic improvement over 
the order 1/5 steps required by the latter. 
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Combinatorial optimization problems (COPs) are im- 
portant in almost every branch of science, from com- 
puter science to statistical physics and computational bi- 
ology [l| . Each instance of a COP requires that we min- 
imize some objective function over a search space con- 
sisting of d configurations. The search space may have 
additional structure, such as that provided by a graph, to 
give a notion of locality. Because d is typically exponen- 
tial in the size of the problem instance, finding a solution 
by exhaustive search is hard in general. One can exploit 
the notion of locality to find solutions more quickly, but 
the presence of many nonoptimal local minima often pre- 
vents efficient convergence to a solution. Therefore, more 
efficient optimization strategies are desirable. 

A well known and often used general strategy for solv- 
ing COPs is simulated annealing (SA) Q. SA imitates 
the process undergone by a metal that is heated to a 
high temperature and then cooled slowly enough for ther- 
mal excitations to prevent it from getting stuck in local 
minima, so that it ends up in one of its lowest-energy 
configurations. In SA, the objective function E of the 
COP plays the role of the energy, so the lowest energy 
configuration is the optimum. The annealing process can 
be simulated with a variety of techniques. Here, we focus 
on discrete Markov chain Monte-Carlo (MCMC) as used, 
for example, in statistical physics 01- MCMC generates 
a stochastic sequence of configurations via a Markov pro- 
cess that, in the case of SA, converges to the Gibbs dis- 
tribution at a low final temperature. More specifically, 
the annealing process is determined by a choice of an an- 
nealing schedule consisting of a finite increasing sequence 
of inverse temperatures (3i < (32 < ... < /3p, and by 
an associated sequence of transition rules {Mi, • • • , Mp} 
consisting of stochastic matrices acting on configurations. 
When the structure of the problem can be exploited by 
a good choice of transition rules, the MCMC algorithm 
can outperform exhaustive search. 



One way to characterize the implementation complex- 
ity of SA based on MCMC is to count the number of 
times that the transition rules must be applied before 
converging to the desired final distribution within an 
acceptable error. For simplicity, we consider regular 
annealing schedules with (3k ~ {k ~ l)A/3 and choose 
A/3 = 0{S/Em), where 5 is the minimum spectral gap 
of the matrices Mfc and Em = maxo- l-E^HI- We assume 
that E has been shifted so that £■ > 0. Let 7 be the spec- 
tral gap of E, defined as the difference between the two 
smallest values in the range of E. By adapting arguments 
from Ref. [1| to the discrete-time setting it can be shown 
that if Pp = PA/3 = 0{-f-^\og{d/€^)), then the proba- 
bility that SA does not return an optimal configuration is 
no greater than e. Thus, for a success probability greater 
than 1 — e, the implementation complexity of SA is given 
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Ideally, AfsA is small compared to the size d of the 
configuration space. Since problem instance sizes are 
typically poly logarithmic in d, AfsA = O {po\ylog{d)) is 
considered efficient. Efficient AfsA is obtained, for exam- 
ple, when computing physical properties of the A'^-spin 
ferromagnetic Ising model in an homogeneous external 
field [a|- However, inefficient AfsA is obtained if the ex- 
ternal field is random [6| , making the problem intractable 
due to gaps 6 that are exponentially small in N. The de- 
pendence of the complexity of MCMC on 5^^ is charac- 
teristic of Markov processes and may be unavoidable [3] . 
Thus, finding new methods with better scaling in 6 is 
very desirable. 

Quantum mechanics provides new resources with 
which to attack optimization problems [8|, |9| . Quantum 
computers (QCs) can theoretically solve some problems, 
including integer number factorization and unstructured 
search, more efficiently than classical computers [10[. 
Still, whether a QC could solve all COPs more efficiently 
than is possible with classical computers is an open ques- 
tion. In this Letter we show that QCs can speed up the 
simulation of classical annealing processes. We present a 
method for transforming instances of MCMC-based SA 



into a quantum simulated annealing (QSA) algorithm 
for which the number of times, MqsAi that the transi- 
tion rules are used is 0{{Em hY'^o^id/ e) log d/{e^/5)), 
a quadratic improvement as a function of 5~^ . This im- 
provement is most significant for hard instances where 
(5 <C 1. The dependence on 1/e can be improved to 
polylog(l/e). QSA is based on ideas and techniques from 
quantum walks [111 ] and the quantum Zeno effect, where 
the latter can be implemented by phase estimation or by 
randomization of an evolution period. 

This paper is organized as follows. First we de- 
scribe a "quantization" of a reversible, ergodic Markov 
chain in terms of a bipartite quantum walk. This is a 
similarit y-tr ansformed version of the quantum walk used 
in Refs. [Il| to obtain quantum speedups in search prob- 
lems. The quantum walk is a unitary operator acting on 
the state space obtained by superposition from the con- 
figurations of the COP. We then explain how to transform 
an instance of SA by adapting the annealing schedule and 
applying the Markov chain quantization. Finally, we an- 
alyze the complexity of QSA to determine the speedup 
over SA. 

Quantum Walks and Markov Chains. Discrete-time 
quantum walks were introduced as_ the quantum ana- 
logues of classical random walks 



1^ 



We focus on the 
bipartite quantum walks defined in Refs. [ll| . 

Consider a d-configuration classical system S with en- 
ergies E\a\ for configurations a. Denote the space of 
ground configurations (minimizers of E) by Sq. Con- 
sider an ergodic, reversible Markov process on S with 
transition probabilities p{a'\a) — viaa' and stationary 
distribution 7r°'. Reversibility is equivalent to the de- 
tailed balance condition 'iT'^m„„i = n'^ rria-'a- Let Ti. be 
the quantum state space spanned by orthonormal states 
\a) for configurations a of S. In SA, tt'^ = e-f^^^^^/Z 
with Z — ^^ e'^^^"^ is the Gibbs distribution at some 
inverse temperature (3. We assume not only that we have 
a classical algorithm to efficiently sample from the dis- 
tribution TOcrcr' giveu a, but also that we have an effi- 
cient quantum algorithm that computes the transforma- 
tion de&red by |cr)|o) i—>- \a) J^a' -\/™ctct^|o'')j with |o) an 
efficiently preparable state of H (e.g. a computational 
basis state). Note that this stronger condition is usu- 
ally satisfied, because for given cr, maa' is non-zero for 
only polynomially many a', and the non-zero m^ra' can 
be computed efficiently by a classical algorithm. 

The bipartite quantum walk is defined on the tensor 
product Ha <8 Hb of two copies of H. Following [11|, we 
define isometrics X and Y that map states of H to states 
of T-Ca «) T-Cb by 



a' 



(1) 
(2) 



nia-a'- From the detailed balance condition, X'^'F = 
Dtt MDtt is symmetric. It follows that X'^Y and 
M have the same eigenvalues Aq = 1 > Ai > . . . > 
A, 



d-l 



> 0. Let \<f>j) be the A^ eigenstate of X^Y. Then 
\(f>o) — ^^ ^/iT^la), which upon measurement in the basis 
\a) has the same probability distribution as the station- 
ary distribution of the Markov process. 
Define unitary operators Ux and Uy by 



UxW)\o)=X\a), UY\o)\a)=Y\a) 



(3) 



Let Dtt be the diagonal matrix with entries tt'^ on the 
diagonal. Let M be the matrix with entries A/o-'o- = 



with arbitrary action on other states. Let Pi and P2 be 
the projectors onto the subspaces spanned by {|cr)|o)}cr 
and {[/jj-C/ylo)!^)}^, respectively. The reflection opera- 
tors through Pi are defined by Ri ~ 2Pi — 1. A step 
of the bipartite quantum walk W based on AI is given 
by W = i?2^i- This walk is related to the one used in 
Ref. [ill by a unitary, but 7r°'-dependent, similarity trans- 
formation, which helps avoid amplitude leakage when W 
changes in QSA. 

The spectrum of W is directly related to the spec- 
trum of M [11|. Define phases fj = arccosAj, so 
that X^Y\(j)j) = cos (pj\(f>j) The spectral gap of M is 
6=1- Xi < (<^i)V2. From Eq. ©, 

PiC/lc/y|o)|(/),)=cos^,|0,)|o) (4) 

P2|0,)|o)= cos <^,C/ic/y|o) !</),), (5) 

so W preserves the (at most) two-dimensional subspace 
spanned by {|0j)|o), C/j(-t/y|o)|(/f)j)}. In terms of the 
Blocli sphere defined by states in this subspace, for j > 1, 
W acts as a Aipj rotation along an axis perpendicular 
to the Bloch-sphere directions spanned by the defining 
states Q^. Thus, the eigenphases of W in this subspace 
are ±2ipj. The eigenphase-0 states are either the quan- 
tum stationary state \4>o) = |(/)o)|o) or orthogonal to both 
Pi. The goal is to prepare \ipo) so that we can sample 
from the stationary distribution of M by measuring the 
ffi'st system. (The preparation of \(f>o) and its relation to 
statistical zero knowledge was studied in |14|.) 

To compare a quantum algorithm based on uses of W 
to the classical Markov chain algorithm, note that W is 
readily implemented in terms of four quantum steps, each 
of whose complexity is closely related to the steps of the 
classical Markov chain, given our assumptions. For the 
purpose of asymptotic comparison, it therefore suffices 
to consider the number of quantum steps W versus the 
number of classical steps based on M. 
Quantum Simulated Annealing. We assume that for 
E^ny /3 > 0, there is a transition matrix M^ satisfying the 
assumptions of the previous section and with stationary 
distribution n"^ = e-^^^^^^Z. Like SA, QSA is based 
on an annealing schedule that we choose to consist of 
equally spaced inverse temperatures Pk = {k — 1)A/? for 
k = 1, . . . ,Q. Let Wk be the quantum walk step opera- 
tor for Mp^, Ii^q) its quantum stationary state {quantum 
Gibbs state for Pk) and (pi,k its phase gap. The goal 
of QSA is to sequentially prepare I'^g'*'"'^) from |-0q) by 



means of an approximate projective measurement onto 



1^1 



fc+i\ 



[l5l | realized by a simulated measurement onto 



the eigenbasis of Wk+i- We assume that the uniform su- 
perposition \^pQ) can be prepared efficiently. If the states 

\^Pq) change slowly enough, the state IV'o ) '-^'^ ^^ °^" 
tained with high probability of success, due to a version 
of the quantum Zeno effect. If /3q is sufficiently large, 
{i/Jq) is a good approximation of a uniform superposition 
of the ground configurations of S, so that we can ob- 
tain such a ground configuration with high probability by 
measurement. The complexity of QSA is dominated by 
the complexity of the simulated measurements, for which 
we give two strategies, one based on the phase estimation 
algorithm (PEA) and the other on randomized applica- 
tions of Wk- Both strategies' complexities are dominated 
by 1/ipi^k- The quadratic quantum speedup is due to the 
quadratic increase of ifi^k over the eigenvalue gap of Mk- 
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FIG. 1: (a) Phase estimation algorithm. The top p-qubit reg- 
ister encodes a p-bit approximation to an eigenphase of Wk^\ 
on readout. The second register's states are in 'Ha®'H,b- The 
first register is initialized with Hadamard gates to an equal 
superposition state in the computational basis. A sequence 
of 2^^ — 1 controlled VKfc+i operations is applied, and the first 
register is measured after an inverse quantum Fourier trans- 
form. If the measurement outcome is |0)®'', the second reg- 
ister is approximately projected onto a 0-phase eigenstate of 
Wfe+i. (b) Randomization procedure. If the PEA's outcome 
is ignored, the overall effect on Ha ® Hb is equivalent to the 
one induced by initializing a set of p bits (first register) in a 
random state r, with r £ {0, ■ • ■ , 2*" — 1}, and by acting on 
Ha ® 'Hb with {Wk+iY ■ Here, double vertical lines indicate 
classical control. 

The use of PEA in QSA is depicted in Fig. ^s.). QSA 
does not need to use the result of the phase estimation, 
though the result could be used to terminate and restart 
the procedure if the measurement outcome is not |0)®p. 
The decoherence it induces in the eigenbasis of Wk+i suf- 
fices to achieve the required Zeno effect. Thus, the effect 
of the PEA on T-La^'Hb is equivalent to the one obtained 
by the action of r Wk+i^, with r chosen uniformly at 
random from to 2^ — 1 (Fig. [IJb)). To exponentially 
reduce the error due to remaining coherences between 
IV'o^^) and orthogonal states, we repeat the random pro- 

cess s times, resulting in a total action of VF^ , '"^ ' with 
< Tg < 2^ — 1 independently random. To prevent ex- 



cessive amplitude leakage into undesirable 0-eigenphase 
eigenstates of Wk, we decohere the second register after 
each randomization step. That is, we measure T-Lb in 
the computational basis and discard the result. The to- 
tal complexity of QSA is given by 0{Q2'Ps) walk steps, 
where Q, p and s are chosen to ensure sufficiently high 
probability of success. 

Let p'^ denote the state after the fc'th randomization 
and decoherence step. We have p^ = |V'o)(V'ol- ^^~ 
sumc that |('(/'o^^|V'o)P > 1 - A^^ for aU k. By ex- 
panding to lowest order in A/3, one can verify that p = 
0{AI3Em)- We show by induction that for 2^ > 2^tt/V2S 
and .s > 1 + log2(2fc)/2 = 0(log(fc)), (V'^|p'=|V'o> > 
1 - 2fc^2 rpi-j^g^ if ^2 < e/(4Q), pQ is the quan- 
tum Gibbs state for f3 = QAf3 with error probability 
at most e/2. We can write p'^ = (1 — x)\i'o){i'o\ + 
i/k{\^o) {ipj_\ + H.c.) + XP±^ where |'0'[) is a unit state 
orthogonal to Ii/jq), P± is a density matrix with support 
orthogonal to |'0o)' ^^'^ X ^ 2fc^^. To make the induc- 
tion argument possible, we add the induction hypothesis 
i/fc < p/2. The induction hypotheses apply to p^ by def- 
inition. Note that (V'S"^^|p*'+i|V'S"^^) : 
We can estimate {4>o^^\p''\iJo ) > (1" 
2^fe|(^o"''l^o')IK^o+'lV'l>l + (^o'+^|p^|^'^+^> > (1 
2kp^){l - p^) - 2vkp > I - 2{k + \)p^. This estab- 
lishes the main induction hypothesis for k + 1. Before 
the randomization step, the density matrix's transition 
between |^o^^) and the orthogonal subspace can be writ- 
ten in the form v'{\tpQ^'^){4>L\ + H.c.) with unit state 
|(/)j^) orthogonal to |i/'o^^) and the other 0-eigenphase 
eigenstates of M^'^"''^, because the decoherence step en- 
sures that the support of Pi is preserved by the oper- 
ator p^ . The estimate on {'^q'^^\p^\'4'q'^^) implies that 
v' < \/2{k -f 1)^ by positivity of p^ [ia|. Because 
|?/;q^^) is stabilized by IFfc+i, the transition is trans- 
formed by randomization to v"{{ijJq^^\\(I)'^) + H.c) with 

v"\<j)'^) = i^' (^ Erlo^ ^+i) I'^-l)- In the eigenbasis 
of Wk+i, the entries of \4>±) are multiplied by terms with 
absolute values 
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< 
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2P-3| 1 ) ^ ^ 7 since the relevant eigenphases 2(p sat- 
isfy 7r/2 > \ip\ > V26. Thus, the choice s = l-|-log2(2(fc-|- 
l))/2 ensures that i'" < p/2. Because the decoherence 
step preserves \ijJq^^), we have Vk+i < i^" < p/2. This 
completes the induction step of the proof. 

To determine the order of the number of quantum 
steps AfqsA required by QSA, let /3/ be the desired fi- 
nal inverse temperature, so that A/3 = Pf IQ- Choose 
Q to be a sufficiently large multiple of (S'iE^/e. For op- 
timization, we let /3/ = ln(d/(2e))/7 = 0{\og{d/e)/-/). 
According to the bounds at the beginning of the pre- 
vious paragraph, this ensures that after measuring the 
final state, the probability of finding a non-optimal con- 
figuration is at most e, with a contribution of e/2 from 
the probability of being orthogonal to |V'o ) and e/2 from 



the Gibbs distribution's probability of not being opti- 
mal. Because 2^ = 0[l/V5) and s = 0(log(Q)), we find 
that Mqsa = 0{Q\og{Q)/sf8) with Q = OiPJEl/e) 
and (}f = C'(log(d/e)/7). ff we anticipate that Q > d, 
we can just search every configuration classically to find 
the optima, so we can bound log(Q) < log{d) to simplify 



Afc 



QSA 




\og\d/e)\ogd 
eVS 



(6) 



The dependence of J^qsa on 1/e can be improved to 
polylog(f/e) by repetition of QSA with an initial tar- 
get error e = f/2 in Eq. ([6]). For optimization, it suf- 
fices to repeat QSA 0(log(e)) many times. Another ap- 
proach that may be used to prepare the desired station- 
ary state with high probability of success is to apply a 
high-confidence version of the PEA [fJl at the end of 
QSA to project onto \iPq), the stationary state for inverse 
temperature f3f. If the projection fails, the algorithm is 
repeated. 

Although the dependence oi Nqsa on Em/ J is worse 
than the one appearing in classical SA, it is worth noting 
that unlike the inverse spectral gap 1/5, in many impor- 
tant applications this parameter is bounded by a constant 
or a polynomial in instance size. 

Conclusions. We presented a quantum algorithm based 
on a "quantization" of simulated annealing algorithms 
implemented with MCMC methods. This quantum sim- 
ulated annealing (QSA) algorithm forces the state to 
closely follow a superposition with amplitudes derived 



from finitc-temperaturc Gibbs distributions. This is 
accomplished by either an explicit measurement using 
phase estimation with quantum walk operators, or by de- 
coherence using random applications of these operators. 
QSA can be used both for combinatorial optimization 
and for sampling from a Gibbs distribution for statistical 
physics applications. In contrast to SA, which scales with 
0(1/(5), where S is the minimal spectral gap of the tran- 
sition matrices, QSA scales with 0{l/\/6). Although in 
general the QSA does not yield a polynomial-resource al- 
gorithm, it reduces required resources by an asymptotic 
exponential factor for the ubiquitous hard cases, where 
the gap becomes exponentially small in the problem size. 
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